{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "from sklearn import mixture"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 导入数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Index(['product_id', 'star_rating', 'review_num', 'review_length',\n",
      "       'review_good_bad', 'svd0', 'svd1', 'svd2'],\n",
      "      dtype='object')\n",
      "(752, 7)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'\\n(product个数,8个维度)\\n0\\tproduct_id\\n1\\t该product_id的所有评论 的 平均star_rating\\n2\\t该product_id的所有评论 的 数量\\n3\\t该product_id所有的review_headline和review_body的长度\\n4\\t该product_id所有的review_headline和review_body的好词个数 / 坏词个数\\n# 4\\t该product_id所有的review_headline和review_body的好词个数\\n# 5\\t该product_id所有的review_headline和review_body的坏词个数\\n5\\t该product_id的所有评论 的 svd0\\n6\\t该product_id的所有评论 的 svd1\\n7\\t该product_id的所有评论 的 svd2\\n# '"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "# market = \"hair_dryer\"\n",
    "# market = \"microwave\"\n",
    "market = \"pacifier\"\n",
    "\n",
    "tttt = pd.read_excel(\"./2c_output/successful_failing_product_2c_\" + market + '.xlsx', 'Sheet1')\n",
    "# print(tttt.head())\n",
    "num_zhibiao = len(tttt.columns)\n",
    "X = tttt.values[:, 1:num_zhibiao+1]\n",
    "# print(X[1:5, 1:7])\n",
    "# print(X[1:5])\n",
    "print(tttt.columns)\n",
    "print(X.shape)\n",
    "\n",
    "'''\n",
    "(product个数,8个维度)\n",
    "0\tproduct_id\n",
    "1\t该product_id的所有评论 的 平均star_rating\n",
    "2\t该product_id的所有评论 的 数量\n",
    "3\t该product_id所有的review_headline和review_body的长度\n",
    "4\t该product_id所有的review_headline和review_body的好词个数 / 坏词个数\n",
    "# 4\t该product_id所有的review_headline和review_body的好词个数\n",
    "# 5\t该product_id所有的review_headline和review_body的坏词个数\n",
    "5\t该product_id的所有评论 的 svd0\n",
    "6\t该product_id的所有评论 的 svd1\n",
    "7\t该product_id的所有评论 的 svd2\n",
    "# '''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GaussianMixture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(752,)\n",
      "(752,)\n",
      "(752, 9)\n",
      "\n",
      "   star_rating  review_num  review_length  review_good_bad      svd0  \\\n",
      "0     4.350632   30.786701     242.221811        12.168524  0.184697   \n",
      "1     4.287915    2.043851     202.662097         1.495319  0.009788   \n",
      "\n",
      "       svd1      svd2  \n",
      "0 -0.122763 -0.019254  \n",
      "1 -0.007730 -0.004008  \n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Gaussian Mixture Model Ellipsoids\n",
    "https://scikit-learn.org/dev/auto_examples/mixture/plot_gmm.html#sphx-glr-auto-examples-mixture-plot-gmm-py\n",
    "\"\"\"\n",
    "# Fit a Gaussian mixture with EM using n_components components\n",
    "n_components = 2\n",
    "gmm = mixture.GaussianMixture(n_components=n_components, covariance_type='full').fit(X)\n",
    "\n",
    "# Estimate model parameters using X and predict the labels for X.\n",
    "# 用X估计模型参数，并预测X的标签。\n",
    "label_X = gmm.fit_predict(X)\n",
    "print(label_X.shape)\n",
    "# print(label_X[1:20])\n",
    "\n",
    "# Compute the weighted log probabilities for each sample.\n",
    "# 计算每个样本的加权log概率。\n",
    "score_X = gmm.score_samples(X)\n",
    "print(score_X.shape)\n",
    "\n",
    "# 输出表\n",
    "X_out = np.concatenate(\n",
    "    [X, label_X.reshape(-1, 1), score_X.reshape(-1, 1)], axis=1)\n",
    "print(X_out.shape)\n",
    "X_out_columns = [tttt.columns[1], tttt.columns[2], \n",
    "                 tttt.columns[3],tttt.columns[4],\n",
    "                 tttt.columns[5],tttt.columns[6], tttt.columns[7],\n",
    "                 'label', 'weighted log probabilities']\n",
    "X_out = pd.DataFrame(X_out, columns=X_out_columns)\n",
    "X_out.to_excel('./2c_output/X_label_score_samples_2c_' +\n",
    "               market + '.xlsx', index=False)\n",
    "\n",
    "\n",
    "# 计算 n_components 个类的中心即均值\n",
    "'''\n",
    "# 这个笨方法也行\n",
    "class_centers = []\n",
    "for ii in range(n_components):\n",
    "    class_centers.append(np.mean(X[np.where(label_X == ii)], axis=0))\n",
    "class_centers = pd.DataFrame(class_centers, columns=tttt.columns[1:9])\n",
    "print(class_centers.head())\n",
    "class_centers.to_excel('./2c_output/n_components_centers_2c_' + market + '.xlsx', index=False)\n",
    "'''\n",
    "print()\n",
    "class_centers = pd.DataFrame(gmm.means_, columns=tttt.columns[1:num_zhibiao])\n",
    "print(class_centers.head())\n",
    "class_centers.to_excel(\n",
    "    './2c_output/n_components_centers_2c_' + market + '.xlsx', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 散点图\n",
    "https://scikit-learn.org/dev/auto_examples/mixture/plot_gmm.html#sphx-glr-auto-examples-mixture-plot-gmm-py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEXCAYAAABCjVgAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfZxVdbn//9c1Mww0iCAMICB3omGEhoKI2o2poY06WorlV8uy85NKK052yhMdK8+PTlq/0tJvoZbH0oMKpY1FRUGYhXEERUUNFRQEUUFgFEduZub6/bHWXqy9Z+89e8/su4H38/Hgway9117rWmtgXftzb+6OiIgIQFW5AxARkcqhpCAiIhElBRERiSgpiIhIRElBREQiSgoiIhJRUpAezcx+b2aXlvH8o8xsp5lVlyuGUjOzi81sUWz7ZDN7LrwP55X7dyLdYxqnIPkws48D/wpMBN4CXgDuAH7i+9E/JjP7FHA78EN3/3Ls9fOA+4A73P1TeR7zv4GN7v6NwkVafma2GGhy9xvLHYt0n0oKkjMzuwq4EfgecCgwFPgscDJQW8bQimUt8DEzq4m99kng2XIEkxJHJRkNPNXdg1Tw9R1QlBQkJ2bWH7gW+Ly7L3D3Nz3wmLtf7O67w/3OMrPHzOwNM3vJzL4VO8YpZrYx5bgvmtnp4c9TzWxF+NlXzewH4et9zOxOM3vdzHaY2SNmNjR8b6mZ/Uv48zgzWxLut9XM7jKzASnn+oqZPWFmzWZ2j5n1yXLZrwBPAmeEnx8InAQ0xY45xszczGrMbKCZbTSzc8L3DjKz583sk2Z2OXAx8NWwmuWBcB83syNix/tvM/t/4/fLzL5mZq8QlFwws7PNbFV4L5aZ2TFZfm9uZl80s3XhPfmemVXleL9GmtmvzWxLuM9N4eufMrO/hT+vBQ4HHgivq3f8dxLuc5mZPWNm283sj2Y2OiW+K8zsOeC5LL8LKRElBcnViUBv4Ded7PcWwbfpAcBZwOfCKpdc3Ajc6O4HA+OAe8PXLwX6AyOBQQSlk7fTfN6A/wKGA+8K9/9Wyj4XAmcCY4FjgE91EtMvwusB+DjB9e9Ot6O7bwMuA241syHAD4FV7v4Ld78FuAu43t0PcvdzOjlvwqHAQIJv45eb2XHAz4GZBPdiLtBkZr2zHOMjwBTgOODcMEbIcr/CNpLfAuuBMcAI4O401zwO2ACcE15X0r0Jf/dfBz4KDAYeAualHOY84ARgQrYbIaWhpCC5qge2untr4oXwW+oOM3vbzN4P4O5L3f1Jd2939ycIHgAfyPEce4EjzKze3Xe6+z9irw8CjnD3Nndf6e5vpH7Y3Z939z+5+2533wL8IM25f+TuL4cP8AeASZ3EdB9wSlhS+iRBksjI3RcB84HFBElxZifH70w78M3wmt4G/h9grrsvD+/FHQRJalqWY1zn7tvcfQNwA3BRGGu2+zWVIFn8m7u/5e673P1vXYh/JvBf7v5M+G/nO8CkeGkhfH9beH1SZkoKkqvXgfp4va+7n+TuA8L3ElUSJ5jZX8Iqh2aCb/X1OZ7jM8A7gX+GVURnh6//EvgjcLeZvWxm15tZr9QPm9kQM7vbzDaZ2RvAnWnO/Urs5xbgoGwBhQ+q3wHfAOrd/e85XMctBA3xt7v76znsn80Wd98V2x4NXBUm4x1mtoPgG/7wLMd4Kfbz+sS+ndyvkcD6+JeALhoN3BiLdRtBCWVEhvikzJQUJFcPE3wjPbeT/f6HoM59pLv3B35K8BCAoGqpLrFjWEUxOLHt7s+5+0XAEOA6YIGZ9XX3ve7+bXefQFCnfzb7qnTi/gtw4JiwCuqS2Lm74xfAVQTJKavwmuaGn/lcvL0gjC1VC7F7QlBdFJf6mZeAOe4+IPanzt1Tq2TiRsZ+HgW8HP6c7X69BIyy7jf+vgTMTIn3He6+LLbPftNrbX+gpCA5cfcdwLeB/2tmF4SNqFVmNgnoG9u1H7DN3XeZ2VTg/8TeexboY0FjdC+Cb99RXbiZXWJmg929HdgRvtxmZh80s6PDB+4bBNVJbWnC7AfsBHaY2Qjg3wpy8fAg8CHgxzns+/Xw78uA7wO/sH1jGF4laJSNWwX8HzOrNrMz6byq7Vbgs2GJzMysb3g/+2X5zL+Z2SFmNhL4EnBP+Hq2+/W/wGbgu+E5+pjZyZ3Els5PgX83s3dD0GHBzGZ04ThSIkoKkjN3vx74MvBV4DWCh9xc4GtA4pvf54FrzexN4Br2NRbj7s3h+7cBmwhKDvHeSGcCT5nZToJG54+HVSeHAgsIEsIzBA/pO9OE+G2CxtRmgiqfX3f7ooO43d0Xh+0QGZnZZIL780l3byMo7ThwdbjLz4AJYVXK/eFrXwLOIUiCFwP3k4W7ryBoV7gJ2A48T+eN5b8BVhIkoN+FcUCW+xXGfw5wBEFD8kbgY52cJ1289xHch7vDKqrVwIfzPY6UjgaviezHzMyBI939+XLHIj2DSgoiIhJRUhARkYiqj0REJKKSgoiIRHr8BFT19fU+ZsyYcochItKjrFy5cqu7D059vccnhTFjxrBixYpyhyEi0qOY2fp0r6v6SEREIkoKIiISUVIQEZGIkoKIiESUFEREJFLUpGBmPzez18xsdey1gWb2JzN7Lvz7kPB1M7MfWbB84RPhClMiIlJCxS4p/DfBzJdxVwOL3f1IgtWpEjNIfhg4MvxzOfCTIscmIpJRU9MarrxyIU1Na8odSkkVNSm4+18JVlqKOxe4I/z5DoL1WROv/yKcpvgfwAAzG1bM+ERE0mlqWsNFF/2Km29+hIsu+tUBlRjK0aYw1N03A4R/DwlfH0HysnwbSV6yL2Jml5vZCjNbsWXLlqIGKyIHnkWL1tLSsheAlpa9LFq0tswRlU4lNTSnWzYx7Wx97n6Lu09x9ymDB3cYpS0i0i3Tp4+jri5YBryurhfTp48rc0SlU45pLl41s2HuvjmsHnotfH0jyWvJHsa+tWRFREqmsXE88+adz6JFa5k+fRyNjePLHVLJlCMpNAGXAt8N//5N7PUrzexu4ASgOVHNJCJSao2N4w+oZJBQ1KRgZvOAU4B6M9sIfJMgGdxrZp8hWPs1sYj3QqCBYM3ZFuDTxYxNREQ6KmpScPeLMrx1Wpp9HbiimPGIiEh2ldTQLCIiZaakICIiESUFERGJKCmIiEhESUFERCJKCiIiElFSEBGRiJKCiIhElBRERCSipCAiIhElBRERiSgpiIhIRElBREQiSgoiIhJRUhARkYiSgoiIRJQUREQkoqQgIiIRJQUREYkoKYiISERJQUREIkoKIiISUVIQEZGIkoKIiESUFEREJKKkICIiESUFERGJKCmIiEhESUFERCJKCiIiElFSEBGRSNmSgpn9q5k9ZWarzWyemfUxs7FmttzMnjOze8ystlzxiYgciMqSFMxsBPBFYIq7TwSqgY8D1wE/dPcjge3AZ8oRn4jIgaqc1Uc1wDvMrAaoAzYDpwILwvfvAM4rU2wiIgeksiQFd98EfB/YQJAMmoGVwA53bw132wiMSPd5M7vczFaY2YotW7aUImQRkQNCuaqPDgHOBcYCw4G+wIfT7OrpPu/ut7j7FHefMnjw4OIFKiJygClX9dHpwAvuvsXd9wK/Bk4CBoTVSQCHAS+XKT4RkQNSuZLCBmCamdWZmQGnAU8DfwEuCPe5FPhNmeITETkglatNYTlBg/KjwJNhHLcAXwO+bGbPA4OAn5UjPhGRA1VN57sUh7t/E/hmysvrgKllCEdERNCIZhERiVFSEBGRiJKCiIhElBRERCSipCAiIhElBRERiSgpiIhIRElBREQiSgoiIhJRUhARkYiSgoiIRJQUREQkoqQgIiKRTmdJNbMnybACGoC7H1PQiEREpGxymTr77PDvK8K/fxn+fTHQUvCIRESkbDpNCu6+HsDMTnb3k2NvXW1mfweuLVZwIiJSWvm0KfQ1s/cmNszsJKBv4UMSEZFyyWfltc8APzez/uH2DuCywockIiLlknNScPeVwHvM7GDA3L25eGGJiEg55LVGs5mdBbwb6GNmALi72hRERPYTObcpmNlPgY8BXwAMmAGMLlJcIiJSBvk0NJ/k7p8Etrv7t4ETgZHFCUtERMohn6Twdvh3i5kNB/YCYwsfkoiIlEs+bQq/NbMBwPeARwlGOd9alKhERKQs8ul99J/hj78ys98CfdQDSURk/5JzUjCzPsDngfcSlBL+ZmY/cfddxQpORERKK5/qo18AbwI/DrcvIpgHaUahgxIRkfLIJymMd/f3xLb/YmaPFzogEREpn3x6Hz1mZtMSG2Z2AvD3wockIiLlks96Cr2AT5rZhnB7NPB0ccMTEZFSymc9hazM7BB3357ricPurbcBEwmSzGXAGuAeYAzwInBhPscUEZHu6bT6yN3XZ/sT23Vxnue+EfiDux8FvAd4BrgaWOzuR4bHuzrPY4qISDcUco1my3nHYKbV9wM/A3D3Pe6+AzgXuCPc7Q7gvALGJyIinShkUsi4jnMahwNbgNvN7DEzu83M+gJD3X0zQPj3kHQfNrPLzWyFma3YsmVLtwMXEZFAIZNCPmqA44CfuPuxwFvkUVXk7re4+xR3nzJ48OBixSgicsApS/URsBHY6O7Lw+0FBEniVTMbBhD+/VoB4xMRkU50mhTMbGC2P7FdT8v1pO7+CvCSmY2PffZpoAm4NHztUuA3uR5TRES6L5cuqSsJ2gsMGAVsD38eAGwgnD7b3bflee4vAHeZWS2wDvg0QZK618w+Ex5bU2iIiJRQp0nB3cdCtPJak7svDLc/DJze1RO7+ypgSpq3ci5xiIhIYeXTpnB8IiEAuPvvgQ8UPiQRESmXfCbE22pm3wDuJKhOugR4vShRiYhIWeRTUrgIGAzcB9xPMIbgomIEJSIi5ZHPymvbgC+Fo5Hb3X1n8cISEZFyyLmkYGZHm9ljwJPAU2a20swmFi80EREptXyqj+YCX3b30e4+GrgKuKU4YYmISEJT0xquvHIhTU1rin6ufJJCX3f/S2LD3ZcCfQsekYiIRJqa1nDRRb/i5psf4aKLflX0xJBPUlhnZv9hZmPCP98AXihWYCIiAosWraWlZS8ALS17WbRobVHPl09SuIyg99GvCXogDSYYhSwiJVbK6gQpr+nTx1FX1wuAurpeTJ8+rqjnM/d8ZryO1kKomN5HU6ZM8RUrVpQ7DJGSSVQntLTspa6uF/PmnU9j4/jOPyg9VlPTGhYtWsv06eMK9rs2s5Xu3mFWCfU+EulhilWdoNJH5WpsHM9NNzWUJPmr95FID1OM6oRSN2ZK5VLvI5EeprFxPPPmnc8VVxxfsKqjUjdmSuXKZ+6jdWb2H8Avw+1LUO8jkbJobBxf0KqE6dPHcfvtq6J2imI3ZkrlyicpXAZ8m6D3kQF/Rb2PRPYLidJHoRszpefJu/dRpVHvIxGR/GXqfZRzScHMHiCYMjuuGVgBzHX3Xd0LUUREyi2vEc3ATuDW8M8bwKvAO8NtERHp4fJpUzjW3d8f237AzP7q7u83s6cKHZiISDbFGNAl+ZUUBpvZqMRG+HN9uLmnoFGJiGShcRXFk09SuAr4m5n9xcyWAg8B/2ZmfYE7ihGciEg6GldRPDknBXdfCBwJzAr/jHf337n7W+5+g5l9qFhBiojElXqSuANJPm0KuPtu4PEMb18H/KnbEYmIdELjKoonr6TQCSvgsUREsir0qG4J5NOm0JmePQpORKSIesostIVMCiIikkZP6i1VyKTwYgGPJSKy3+hJvaXyWWTnITObY2Znmlm/1Pfd/aOFDU1EZP/Qk3pL5dPQfCnwXuB84Htmtht4yN3/tSiRiYjsJ3pSb6mck4K7rzOztwlGL+8BPgi8q1iBiYjsT3pKb6l8qo/WAvcDQ4GfARPd/cxiBSYiIqWXT0Pzj4ANwEXAF4FLzaxbFWNmVm1mj5nZb8PtsWa23MyeM7N7zKy2O8cXEZH85DPNxY3uPgM4HVgJfAt4tpvn/xLwTGz7OuCH7n4ksB34TDePLyIVpqf0189VuuvpydeY88prZvb/ETQ0HwQ8TDAh3kPuvq5LJzY7jGAivTnAl4FzgC3Aoe7eamYnAt9y9zOyHUcrr4n0HIn++om1oOfNO79H1LNnku56gB5xjZlWXsun+ugfQKO7v9vd/8Xd7+hqQgjdAHwVaA+3BwE73L013N4IjEj3QTO73MxWmNmKLVu2dCMEESmlntBfP59v+emupydcYzb5JIVfAR8ys/+AYD0FM5valZOa2dnAa+6+Mv5yml3TFmPc/RZ3n+LuUwYPHtyVEESkDCq9v36+I4/TXU+lX2Nn8hmncDPBt/pTgf8E3iRIFMd34bwnA41m1gD0AQ4mKDkMMLOasLRwGPByF44tIhWq0vvrp/uWny3GTNdTydfYmXySwgnufpyZPQbg7tu72jvI3f8d+HcAMzsF+Iq7X2xm84ELgLsJBsv9pivHF5HKVcn99fv375N1O51011PJ19iZfKqP9ppZNWGVjpkNZl97QKF8DfiymT1P0MbwswIfX0Qko+bmXVm3e3KvolzlO07hPmCImc0B/gZ8p7sBuPtSdz87/Hmdu0919yPcfUa4qI+ISElkaw/oSTOddkc+01zcZWYrgdMIGoXPc/dnOvmYiEiPka3NI9/2hp6q06RgZge7+xtmNhB4DZgXe2+gu28rZoAiIqWUqT1g+vRx3H77qmj8QU/rVZSrXEoK/wOcTTCKOd5F1MLtw4sQl4hIUTQ1rUlbEsj0ekK2UkRnn+1Jch7RXKk0ollEcpVpRHV3Rlr31FHa3R7RbGa/MbOLzKyusKGJiJRGptHGnY1CztbrqKePYE6VT++jHwDvA54xs/lmdoGZdd6JV0SkQmTqXdSdXkc9fQRzqnxmSX3Q3T9P0IZwC3AhQcOziBwgytlPvxDnTrQLXHHF8UnVPI2N45k1axoTJw5h1qxpnfY6SnfMhoYjOOWU0V2OrVLkM6IZM3sHwWymHwOOI5jlVEQOAPG689tvX1XSuvNCnjtd76KmpjXccMM/aGnZy7p12znhhBHRPrn2Olq6dD0tLXtZunR9j2lXSCefNoV7CNY+OJVgHqRx7v6FYgUmIpWlnHXnxT53tuNnKl2UMr5SyqdN4XaCRPBZd1/i7oWe4kJEKlg5686Lfe7Ojt/YOJ6bbmrI+O1/f2pXyGeRnTqCxXBGufvlZnYkMN7df1vMADujLqkipVPO/vjFPnd3j9/Txipk6pKaT1K4h2AA2yfdfWLYvvCwu08qbKj5UVIQEclfIVZeG+fu1wN7Adz9bdIvjCMiIj1UPklhT1g6SEydPQ7QLKYiIvuRnLqkmpkBPwX+AIw0s7sIVk/7VPFCExGRUsspKbi7m9mXgOnANIJqoy+5+9ZiBici+5dSNMb2tAbfSpPP4LV/AIe7+++KFYyI7L9KMfitnAPs9hf5tCl8EHjYzNaa2RNm9qSZPVGswERk/5JtgFehps8oxCCyYk3lUajjFnuqkXySwoeBcQQjms8hWGPhnGIEJSKVo1APoUwDvAq5zGV3B5EVa8nNQh23FEuC5jMh3vp0fwoekYhUjEI+hDJNF1HIKSJymZIim2JNV1GoUlIpptPIp6QgIgeYQj+E0k0XUegpIjqbkiKbrsSSy0O9UKWkUkynkdcsqSLSs+XbM6cU6xJnW+ay1PKJpalpDXPnruDPf36BPXvasjZsZzpuuqSb7ZyluFdajlPkANHVZSPVxbOj+L2Mu+KK47nppoYuHafUS3lmmuZCJQWRIqjEB2m+30oT0q0/UGiVeL+yid/LhK6UpCqplJSgNgWRAitFD5GuqNTpnSv1fmUzffo4eveuBqCmpoqGhiO6/C2/O20gxaCkIFJglbrgSnd75hRLpd4vyN6InKh5r6oy6upqmT17CbNnL8n585VK1UciBVaKxtmuKkVVUL4Kfb8KVRWVbXT0okVr2bOnDYA9e9pYsOBpAFavDpatnzPn1B47ulolBZECq9Rv5JWqkPerkFVR2Uow8ao4S1lAIHHOSi4BZaOkIFIElVZPXOkKdb8K+SDO1gYTT2Tnnz8h6XOJa6jUNpzOqPpIRMouW5VPPtVBnVVF5XOedD2D0u0DMHv2Epqa1tDYOJ45c04FKrNnUS7KMk7BzEYCvwAOBdqBW9z9RjMbCNwDjAFeBC509+3ZjqVxCiI9W7a++l3px5/pwZ/reXr3rua008Yyc+aUnD/fExViOc5CagWucvd3EazPcIWZTQCuBha7+5HA4nBbRPZj2ap8ulIdlKkqKtfz7N7dxsKFz3dok+ipbQT5KktScPfN7v5o+PObwDPACOBc4I5wtzuA88oRn4gUVraumdnq3gtZL5/reRKyNS7nOy/S7NlLekzX1LJPc2FmY4C/AhOBDe4+IPbednc/JM1nLgcuBxg1atTk9es1WatIpcql2qVQbQq5xJLtPPG5jNLFmvh8//59aG7elbaKKnF8oMNUGIWodirU/chUfVTWpGBmBwEPAnPc/ddmtiOXpBCnNgWRzCph+ogrr1zIzTc/Em3nOz9QqXV2zzIludTXTzllNAsXPt/h8925/kK2a1RamwJm1gv4FXCXu/86fPlVMxsWvj8MeK1c8Yn0dJUyfURP65rZWffYTG0Lqa8DHaqksl1/LqOfS9GuUZYuqWZmwM+AZ9z9B7G3moBLge+Gf/+mDOGJ7Be6OgFeofXUrpmZxLu9AvzhD8/T1LSmQ3fYmTOnMHPmlKzVTQm5jn7u379P1u1CKFdJ4WTgE8CpZrYq/NNAkAw+ZGbPAR8Kt0WkCyrpG/r+NJivsXE8DQ1HRttr127n/PPvBWDWrGlMnDiEhoYjo2/xN93UwJw5p3ap9JFaemhu3pX0udTtQih7Q3N3qU1B9leFaA+ohDaFfCQae4EO4wQqQeJ+3nvvU2zZ0pL0XkPDESxduj6pYbl372quuuqkrKWExHFT2wqAnF4rdJuCkoJIBdrfBkolpBv5m9DUtIYZM+ZHE8317l3NvffO6DCSGChooss1ccZ/J2b7ZkkFqKqCM888Im3DcnW10dbm1NRUMXx4Py655JgO154ujkwN9NnuYT60yI5IBensQdSd9oBKLR3Mnr2E73znISB5NtGE+MyjEAwiS1SjJB7Gt976KGbBe6l171257lzr8pua1jB79pLod+IeJIL29mBCvKuvfh8nnDCiQ0kBoK0tyB6tre1s2NAc3YPEtcfjjvdKSjdlR1PTGm644R+0tOxl3brtnHDCiIL/jjUhnkiJ5dIrqKvtAZXS4yhVU9MafvrTRzq8Fjd9+jhqa6uj7d69q5k+fVxSgtyzp43du4PEkVr3nrjuCy+cz1ln3ZXTtedSl584diKRQfA7+ehHJzB6dH+OO25Y9HCeN+98GhqOoKYmeLTW1FRFi/HE3XXXE9GgtkTcM2Ykx51u9tj9tveRyIEsl1JAV3vsVEqPo7hM6xmnu+b582ekbVNIfGOura2OSgrxZJlumoqlS9d3Wu2W6dt4vPRwyimjk2I3g0mTDuW++56hrc1Zv76ZGTPmM3/+jGiSvNTqrmuvfZCVKzdHx9i06U1uvvmRqGoJgoSXGnfq+helWKtDSUGkxHL9j92VBXG6+9BIVwXT3eqo1PWMBw7sw2c/e3za+vB015yaIBPHjM9c+sIL2+nduzoqRUByUsx0DemS75VXLkxKrKtXb4mqiiCoOlq27KWkGPfsaUtKwKnX0dg4PmoL2Lr1LV555S2AqK2htbU9bdzp7s+sWdOiNoViJHw1NIuUQTHr/bt67Fx7wOQbb3cazfMZXVxbW83RRw/hySdfS5qmIt9rSFeyqaoy2tszPytra6ujkkK2uGfMmB+t0gZB9dJXv3oyq1ZtZvHiF6ISULa2jWL3PlJJQaQMirksZlePnam+urvVUV2tCsulETi1vWHatMO45poPZP3m39k1JOKdPXtJ1I6QKSFUVcGQIX3p3buG5cs3dZjuIh53U9OapIQAMGJEP044YQTNzbu46qphHbqupiaXUlQPKimICJC56qkQddhdSVS5PAD79+8T1csn4uusHr5//z6cddZdQHK7RWo7wKhRB/Pss6+zZ08btbXVtLa2JyUHM+jduyaqCkr0Kmpu3pU27nSNwscfPyKnNR4SyaUUI5qVFET2I7lUt2R6P9M3+ngdNgQT3HXl+PnK1j6SGOS2ePELUb18fBRx/NzLl2+ivr6OwYPrOOOMI/j+95dFXV8XL36Be++dAaTv9tq7dzUNDUcwadIwrrvub0nxucPbb7d2uP45c05NG3fq9BgXXDCBoUP7Zkx8ufQ00ojmNNSmIJWoHGMFOqtv7uoqZvEVydxJO610LtNOpx431wFjqY3K8fPExUsMiXPHx0YATJ48LKkXEASDwoCkgWJxEycOYdSog9MOTEt10kkjOfbYQ3OaWjvdzKqZSgq1tdWcfvpYJk0aFo1TUJuCSA+R64CoQuusuqUr9dGpXT0TWlr2Mnfuig4Ptvj7mY6fy/1JN6ArU9dWIKkHT2qvo7jnntuW1GicGAsBcOutj3ZINBAMtHv66S0degmls2zZSyxb9lJSI3dqVVUiriuvXEj//n045ZTRSe8nJEpuiRJRorvqrFnTOp02ozuUFKQsyj3qtpjnL9dYgXj1RG1tNS+8sD2p2ieX7qqp9yW1zjr+QF28+IVo/9QHdbb2h9T7kxiXkDjv8uWbuO66v9HW5klJY+7cFWkTwuTJwzjjjCOSvkH379+HK69cyMEH907a9403dofXAcceO4xrrvlAdH9OP31sUmlg4MA+bNsWVM8E1+xMnjyMxx9/tdPkkLiueIkmXVVV/H7NnNnhSzsAGza8kTRgr7l5V1HXo9CIZslZLvO953qcco66Lfb5yzU7aXxErRkd1hlON0I2Lt19Sa2zHjny4OjnxDQU8etN1MFnKx2ljlxetGgdM2bMj0b1JhICBA/BL3zh90yZcguLFq1Le7xp0w5jzpxTo2traDiS66//Ozff/EiH8QQJ7e3B5+JVNUA0Erm2tprPfvb4aDvxmaFD+3LIIcmJ0oyk64F96yikTtsxe/aStMktXZvB7NlL+OhH7+kwkiOpxOsAABggSURBVDqR8Ir1/0YlBclJIatEyj3qttjn72oXzEKde9GitezeHXzjTb2+bL2A0t2X1JLC8cePYMuWlqTSRvx6E3Xpy5dvSnv9iZLF0UcPier249+601XfbNjQzIYNzUmvJQaTxUchJ87/k588Eg00y6S2dl+1UbpqKTM44YQRfPWrJ0dJqqoK/vCHtR16IJ1//gTe+c5BNDWt4aij6hk6tG907NS2j9WrX+PZZ1/vMNAuXUP69df/PUqOELRtNDaOj0pExaqaVFKQnBTyQVqKofrlPn8xxyF0Jp/ri1cXpftc6rfXoUP7pk14ib9TH66p/fTjDaeJB2Pv3tW0tTmtre3U1FRRXW3s3t3WYSbShNraakaOPJi9e9u55JJjAKLZVTN9Jj6dBARVRel6+STs3t3Gtdc+yIoVlwNw/fV/D5NX8sHd4YEH1kQN8M88s4VJk/Y1NH/lKyexatVmnnpqC+vXB4ltz542GhqOYOzYQzI2SC9atDYpWVZXG3PmnKpxClI5CvkgLec36Uo4f7Hlen3pSn/pPpf6e8+U8NI9XFP76ccHmsUfjNdf//foM4n1BxKvJx6ONTVVvOc9Q3n88VdZu3Y7AN///jKOPnpI9G08XUIYPbo/xx8/ImngWF1dbdS1dvr0cdx226NJ39wBVq7czOzZS2hu3pW1DSH+ubY2Z+XKzVEpKN7gnEhcVVUwadKwrNNex/+/JUY9Z/p9FJq6pErOyt04LIWVab7+VF1ZbyAh28L2idfPOuuupAbehoYj+N3vLqapaQ0XXjif3buDB+nVV7+P5uZdHbqODh5c12HBm7hE8rn5sXUweQCs3AHLXo/imzVrWlLyiRs1qj8//vGHozi64oorjmf69HF85CP3RFVPNR8YzHsuexdDX9rDzKPH5DXmo1D/D9UlVbqtnFUi0nWZHiKFnJgvcY6GhiOj2UNraqqYNWtaUhVTupLIq6++lXSsxHbQNhI8iNvbiUoPcVVVxvbtu5K243X+VVXGzJlT+OW6TXD2u6BPNTQMg/98Bpa9TkvLXn70o+UZSwKDB9fR2DieiROHdBjfkIvEfZ07d8W+uE4aROvV41nZx2BIDX++fgnzST9rbLr7Xuz/h+p9JFIETVu3cuWzz9K0dWv2/brZo6uz82TradVZb6ScY4id4/77/5m0qEwuI26HDu0Ll42BWyfDZWOCbZJ7KSUahletSn4wDxz4jqQH+plnjuPrX38fo0f3Z/LkYdx338dobBzPPw9qCxICBH9PHhB9ZufOPWnjqqmp4pprPrAvxpjqast6TWYk9cJKSnyTByTFsmdiv7zWRZixejWDHnqIGatX5/yZfKikIFJgTVu3ctHTT9PS3s7tr7zCvAkTaKyv77hfN3t05XKezhomC/GtM36O1tb2DnMRRfFmuN4d5w+DsRY8Scf1ZccLyb174n+nGj26Py0te6PSTmIAWGp9/VE7q1m9px1qq2BPe1CFlEbqmsoQVLNNmjSMRYvWRQko3mgN0KdPDbt27ZvyIlErv2jRWpYv35TUrZRHd1B9zgjaehnsaqN29ZtMn3V8hrubbMbq1SzYshUMFmzZyozVq5k/cWJOn82VkoJIgS3ato2WsE9kS3s7i7ZtS5sUutuTJJfzxBtR46N3c5VL/XVqo+h55x2V1C0z0aCb6XqfOKgVrDY4mBlPHBTsE68+SoyJmDlzSjTFdO/e1VxzzQdYXrOLppdfo3H4EBobOo69mDt3BYveaoYJ4XsWtBVMHT6UhQufi6bwOO20sR0myDv//HtpbW2nqip7yWDv3rakEc81NVVhnM936Pk0eVctQ/93D6+OrGXoS61MOnHfZHmd/f7/8NrWffU7Fm4XmJKCSIFNHziQ2195hZb2duqqqpg+cGD6/TLU6ee6SH2u50l8a823T0muJZnGxmDhl+uu+xutre0sXPhch3UMbr99FbNmTaOurleH6x26dhc7B/ci0Z906NqgyinTDKj33jtj3/05aRA3PP00LUf04jnfzqo5f4wabmcvfJzrHnyath3bgyqbxACzXlVM/NR45p8xNVr4Jl3p4tprH4we8tnWUoCg5NDQkJxwE43n8YV0eveu5sknX2Plys3U1fXijFnT8hp3sHdnK/Srju7V3p1da/zORklBpMAa6+tpGDiQJdu3c+ohh6QtJUD6htf4g/i22x6N+r+ne2A01tcz67DDaNq6lcb6+oylkUR3zdTVwTqTT0mmuXlX0ijkdGsxNDfvStvQ/LGdB/GdOzfAiYPg4df52KEjowXqMzVYRzO2PvtsVFrabbDw1S0s/c5KZs0/i+trt9F2zjD40BCYvzGoNkpUH63YTtPuNdEDed267dE6ywkbNjTDSYM69FgC0o6HqKur5ROfOIZFi9Z2aDyfOnUEAwb0ZvXqLdFxWyYP4M51m/IqLbbVWFKdWltN9hJMVygpyAGvaetWFm3bxvSBAzM+wPMxe906FoQNvwu2bmX2unXMOfzwtPum1ulnm4Au9YHRtHUrN2zcSEt7O8+t38Cqec8x8+gxcNKg6Hq6M74kn8/muhZDujaMVas2w8IX4ecvBtsNNUlrEmRrsI6XltjVBit30NKyl6aXX6P1iGCqCfpUw5EHsW/gmTPp2ENZtLDjHEzxhLV36gC48vCkHkujNrYyceJggA6zpi5Y8DT33//PqF0FiJLKw4/uoHrRpqDkcdIg+I+gJ9SmvU7Vpjdof2hrTr+fPr2r2ekeZaU+vauz7t8VSgpyQMu1UTjfY6ZuZ0oKqeIP19SpqlMfGPE2hcS35D8vXof1exe7DW7dvJnTxxzCrPln0bxwU9792vMZ5Jdp35wHCca/kZNHd9n6euZNmMC1C59k5dxnom/zR+2sZl1VFS3t7dS0OqPHDWRtbfigrq1mle+mLuXbfKIhOVEqa52U3EuIKQP48deOp7FxPDNmzE8bT+vUQ2DyANpW7gADvhE8/L1hGK1hN9h476O2XgaT+lO97PWk0hCkb89pwZNKCi0UfpyZkoIc0HJtFM5HY309qzdsSNrO+bON2Repj0v3LXnP5AHBwwjY487CbdtYelAV8741tUvX1dg4Pip5EFZTZds31772cXUfGgYTh0FNFTQOp2519oSU+rBsrK9n0f++zcpY9c7Q53cx75KpUYlp7tyVrO1XEzyMd7Xx6sKXWHXfM0lxpE69feYRg1iwpy1oi9jTxkm1/aIqvvtSPgsEie2bE4IqqrOH03fNTt6KJ5VPjIJPj4GW1qDuKVEH9VYbbW2eVBrK1J7TTjuwr3QQbBeWkoIc0HJtrM1HolSQqOvPtZSQkPogzfRQTXxLnrtyHX++fiV7lr1Obe9q7CNBySGhO8kun5JUV6vhHhzlUBM+6GqqeHBUUG2WLqFkelimK1kktbMcPYZF311M6zEHU/PEG7CzJqlHUHW1UV1dlVQqW16zC0gkGuOUD44BgiQd/2xVFRx2WH92XDyGN2rDrkG1VbxVS5Cs+1RjbY6P77cvEcT72vat7lAaytieU5UytCx1uwCUFOSAlniwFrJNAYLEkG8y6IrG+noaz6inaXd/Fk0ISxdHD2Luyy+zePt2drt3K9nlWpLqTjVcc7+qrNtJ8WR4WOZS1VX18DZ4cAvUVPF44sWTBsGUAXxk7HA+cfiIpM/Pvu1PkGiXqK3irtUv0bxwE/3794l6UVVXG1/72nuZM+dUpvxpGSuJDYTbtgd+sQEmD6Dv9OHsPCi8LjNod6gyqve2c8bQwcycd3xOo817tRt7q/aVMnq1q6FZpOAy9dzpSTqULurrC9KAnmtJqjvVcFVtDr0seTtTPFnaGjqbFjzRCysaAR1r8F1YVcUnJgzipvg9HD6E1Xu2RtVHm36zgZv/uiWaLyl1dtMzxh3KyvUbguo7d2rWv03rstepW/UGZ541ngWEbRju8OAWqne28bXTJjBn9ns6xJspyVXtbYc++5JL1V5VH4l0W6F7G1WqQiS7XEtS/Wtqsm5nc1z1O1jmu6Jvv8dVvyNzPF2c4TZ1VToz2B1r8E2XyE6YOoKaJ7fTClRVJS/3mW71s+bW1qg9BzOmf2Q8Y2sGBYlr8iDuX706OBbGmbX9mfmBMR0G26Vea+r1Wa+qpPYI66Xqo4KypUujn/2UU8oWRyGU4kFXu3Qpe4FewJ4Kul+z163Luf6+aetWZjz1FHvcuXXzZua/+91J96tQ97Gr/7bynQEzl2uPH/OXh++Nxk/Ep0fIdt2ff+gJNvWF+5/aSONH9g3wip+7ubU16TOp24k4EktvxkcOHzt6EMtefjnYyYxjRw/Kek9SH5Yz/vgIS95+g6O2Gcc+uotXzxnMI9W7qH9uF9eMGAXA3LkraPvZZBhcy0FvO9N+splX6cNjDu0WlE4e+OZyXm19nJaWsAroi0fSGq7m2V5jVE8bRNuyYIGcF17Yzow/PsI/D3GO2m60/OhZXh3XB84dCNUGDjMnH07jGVOBYDxF4o60G6w5rT8Xt2xm6C9f4mP/rMp5zeXd1RBlHrNwu7AqbupsMzsTuJGgif02d/9utv27OnV2/D9tQk9NDPH63LqqqoJ0q0yVSAgJlZIYZq9bx3diPX2+PmpU1sRw1hNPsHDbtmi7YeBAfndMsFBLoe5jV/9tZZpaOpNcrj1+zKpvT6D9fYOjZ8oF9fXMnzgx63Ufdt8SNg2w6JvpiB3Oxo+c2uHcF9TXszCsQkp37+LTYEMwwd38+TNobByf9X51dk9m/PERFtTu3NeA+2QzHN0/2q66+yWqfr6e1luPhdF99+23/i1q/rGd1o+NjKp7uHNDNFYCgD+/N2jITXymvZ2GH7wSrKZ2yUi4ZNS+9+7cAOccCv1ro9cOajfePC24jnctX84/3347doHsO+/SLXDtMzn9zm3JkuTG5fZ2/NTM6zJkk2nq7IqaJdXMqoGbgQ8DE4CLzGxCeaOqfOnqcwstdbn0jsunl0e6MQFdVYr7mPX8aRpRs8nl2uPHbD+m/77qDWDJ9mChmmzXvakvST1lNvVNf65/trQwb8IErhg+POPEfPHBeInR1Z3p7J4sefuN5J487+yXtN1+/MCg2uewuuT9DqujdcohSdU9nJhcQokSQuL9qirGjj0kaJs4cVDyeycO2pcQwtd2Vu37wv1CPCFA8nmPHZDx+jpInRkw00yB3VBRSQGYCjzv7uvcfQ9wN3BumWOqeNMHDqQu/PZQqG6VqXp1sl0uqQ+fzr7Zzxw+nN7hf6TeZswcPjx6rxT3MZvp08dFC77nMro1l2uPH7Pqieak1SRPPeSQYJ8s1z3iLZImTxrxVuZzN9bXc9M735kxjt6x0bfxNZKz6eyenPqOg5Mnd3r2zaTtqke2UVNTBRtbkvfb2ELNiu2xQc4OD7+edGza25M/0+774nn49eT3Hn4dmvckvXZQrGfQ5H79ko8dP+9jOzJeXwcdYip8Q3NFVR+Z2QXAme7+L+H2J4AT3P3KlP0uBy4HGDVq1OT169d37XxqU8jL/tCmANnvldoUOl73YfctYVPfIEFszNCmkNN9z9CmANnvV2f3JNc2hQc/PYSWQb0Y3l7FefftDMYiTOxF09atHLx6J2/8cA1HHVUftSnMnDmFc/tugqqgC6mfdmpSPK+eMzipTQFgyReHsavWkqqOEk5euZKVb77J5H79GN67N394bStDN7Xm1aYAYIsXB6WY9nb8tNM63T/jcTJUH1VaUpgBnJGSFKa6+xcyfUbLcYqI5K9HtCkAG4GRse3DgJfLFIuIyAGn0pLCI8CRZjbWzGqBjwNNZY5JROSAUVHjFNy91cyuBP5I0CX15+7+VJnDEhE5YFRUUgBw94XAwnLHISJyIKq06iMRESkjJQUREYkoKYiISERJQUREIhU1eK0rzGwL0LUhzYF6oOsT5hSP4spPJcZViTGB4spXJcZViJhGu/vg1Bd7fFLoLjNbkW5UX7kprvxUYlyVGBMornxVYlzFjEnVRyIiElFSEBGRiJIC3FLuADJQXPmpxLgqMSZQXPmqxLiKFtMB36YgIiL7qKQgIiIRJQUREYkcMEnBzH5uZq+Z2eoM75uZ/cjMnjezJ8zsuAqI6RQzazazVeGfa4odU3jekWb2FzN7xsyeMrMvpdmnpPcrx5hKfr/MrI+Z/a+ZPR7G9e00+/Q2s3vCe7XczMZUSFyfMrMtsfv1L8WOKzxvtZk9Zma/TfNeye9VjnGV6169aGZPhufssJpYUf4fuvsB8Qd4P3AcsDrD+w3A7wmW1J4GLK+AmE4BfluGezUMOC78uR/wLDChnPcrx5hKfr/C6z8o/LkXsByYlrLP54Gfhj9/HLinQuL6FHBTGf59fRn4n3S/q3LcqxzjKte9ehGoz/J+wf8fHjAlBXf/K7Atyy7nAr/wwD+AAWY2rMwxlYW7b3b3R8Of3wSeAUak7FbS+5VjTCUXXv/OcLNX+Ce198a5wB3hzwuA08zMKKIc4yo5MzsMOAu4LcMuJb9XOcZVqQr+//CASQo5GAG8FNveSAU8dIATwyqA35vZu0t98rD4fizBN824st2vLDFBGe5XWO2wCngN+JO7Z7xX7t4KNAODKiAugPPDaocFZjYyzfuFdgPwVaA9w/tluVc5xAWlv1cQJPJFZrbSzC5P837B/x8qKeyT7ttIub9ZPUowP8l7gB8D95fy5GZ2EPArYJa7v5H6dpqPFP1+dRJTWe6Xu7e5+ySCNcWnmtnElF3Kcq9yiOsBYIy7HwP8mX3f0IvCzM4GXnP3ldl2S/NaUe9VjnGV9F7FnOzuxwEfBq4ws/envF/w+6WksM9GIJ79DwNeLlMsALj7G4kqAA9WpOtlZvWlOLeZ9SJ4+N7l7r9Os0vJ71dnMZXzfoXn3AEsBc5MeSu6V2ZWA/SnhNWGmeJy99fdfXe4eSswucihnAw0mtmLwN3AqWZ2Z8o+5bhXncZVhnuVOO/L4d+vAfcBU1N2Kfj/QyWFfZqAT4at+dOAZnffXM6AzOzQRH2qmU0l+H29XoLzGvAz4Bl3/0GG3Up6v3KJqRz3y8wGm9mA8Od3AKcD/0zZrQm4NPz5AmCJh62E5Ywrpe65kaCdpmjc/d/d/TB3H0PQiLzE3S9J2a3k9yqXuEp9r8Jz9jWzfomfgelAak/Fgv8/rLg1movFzOYR9E6pN7ONwDcJGt9w958SrAvdADwPtACfroCYLgA+Z2atwNvAx4v9HyR0MvAJ4MmwThrg68CoWGylvl+5xFSO+zUMuMPMqgmS0L3u/lszuxZY4e5NBMnsl2b2PMG33o8XOaZc4/qimTUCrWFcnypBXB1UwL3KJa5y3KuhwH3h95wa4H/c/Q9m9lko3v9DTXMhIiIRVR+JiEhESUFERCJKCiIiElFSEBGRiJKCiIhElBRERCSipCACmNksM6sr8TlPMbOTYtufNbNPljIGkVQapyBCMG89MMXdt+bxmWp3b+tkn5pwYrd0730L2Onu388nVpFiUlKQA044ZcC9BPPEVAPzgdnAGmCru3/QzH4CHA+8A1jg7t8MP/si8HOCKQducve70xx/KbCMYBR2E8HaD98Aagmm3bg4PO4/gDZgC/AF4DTCJBEeYznwQWAA8Bl3fygszfw3cBTBVAtjgCvcvcMCLCJdccBMcyEScybwsrufBWBm/QmmB/hgrKQw2923hdNELDazY9z9ifC9Xe7+3k7OMcDdPxAe/xCCBW7cghW7vuruV5nZT4mVFMzstJRj1Lj7VDNrIJgC5XSCRWi2u/sx4aynqxApILUpyIHoSeB0M7vOzN7n7s1p9rnQzB4FHgPeDUyIvXdPDueI73MY8EczexL4t/B4uUjMBLuSoEQA8F6CmTxx99XAEx0/JtJ1SgpywHH3ZwmmPn4S+C9LWcvZzMYCXwFOC+fP/x3QJ7bLWzmcJr7Pjwmqmo4GZqYcK5vEVM1t7CvVF30VMjmwKSnIAcfMhgMt7n4n8H2CdbLfJFj7GeBggod6s5kNJVjgpDv6A5vCny+NvR4/Z67+BlwIYGYTgKO7GZtIErUpyIHoaOB7ZtYO7AU+B5wI/N7MNocNzY8BTwHrgL9383zfAuab2SaCxuWx4esPAAvM7FyChuZc/F+CKbGfIKjaeoJgyUqRglDvI5EeJGz47uXuu8xsHLAYeKe77ylzaLKfUElBpGepA/4SLk1qwOeUEKSQVFIQ6SIzu5lgLELcje5+ezniESkEJQUREYmo95GIiESUFEREJKKkICIiESUFERGJ/P/xijfd65lvWgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy import linalg\n",
    "color_iter = itertools.cycle(['navy', 'c', 'cornflowerblue', 'gold', 'darkorange'])\n",
    "def plot_results(X, Y_, means, covariances, title, X_pltdimen1, X_pltdimen2):\n",
    "#     splot = plt.subplot(1, 1, 1)\n",
    "    for i, (mean, covar, color) in enumerate(zip(\n",
    "            means, covariances, color_iter)):\n",
    "        # 由于DP不会使用它所访问的每一个组件，除非它需要它，所以我们不应该绘制冗余组件。\n",
    "        if not np.any(Y_ == i):\n",
    "            continue\n",
    "        plt.scatter(X[Y_ == i, X_pltdimen1],\n",
    "                    X[Y_ == i, X_pltdimen2],\n",
    "                    10, color=color)\n",
    "        \n",
    "#         v, w = linalg.eigh(covar)\n",
    "#         v = 2. * np.sqrt(2.) * np.sqrt(v)\n",
    "#         u = w[0] / linalg.norm(w[0])\n",
    "#         # Plot an ellipse to show the Gaussian component\n",
    "#         angle = np.arctan(u[1] / u[0])\n",
    "#         angle = 180. * angle / np.pi  # convert to degrees\n",
    "#         ell = mpl.patches.Ellipse(mean, v[0], v[1], 180. + angle, color=color)\n",
    "#         ell.set_clip_box(splot.bbox)\n",
    "#         ell.set_alpha(0.5)\n",
    "#         splot.add_artist(ell)\n",
    "        \n",
    "#     plt.xlim(-1, 200)\n",
    "#     plt.ylim(-1, 50)\n",
    "    # plt.xticks(())\n",
    "    # plt.yticks(())\n",
    "    plt.xlabel(tttt.columns[X_pltdimen1+1])\n",
    "    plt.ylabel(tttt.columns[X_pltdimen2+1])\n",
    "    plt.title(title)\n",
    "\n",
    "# 绘制好词坏词\n",
    "X_pltdimen1 = 0\n",
    "X_pltdimen2 = 3\n",
    "\n",
    "plt.figure()\n",
    "plot_results(X, gmm.predict(X), gmm.means_, gmm.covariances_,\n",
    "             'Gaussian Mixture ' + market, X_pltdimen1, X_pltdimen2)\n",
    "\n",
    "# # Fit a Dirichlet process Gaussian mixture using five components\n",
    "# dpgmm = mixture.BayesianGaussianMixture(n_components=5, covariance_type='full').fit(X)\n",
    "# plt.figure()\n",
    "# plot_results(X, dpgmm.predict(X), dpgmm.means_, dpgmm.covariances_,\n",
    "#              'Bayesian Gaussian Mixture with a Dirichlet process prior', pltdimen1, pltdimen2)\n",
    "\n",
    "plt.savefig(\"./2c_output/scatter \" + market + \".png\", dpi=500,\n",
    "            bbox_inches = 'tight', pad_inches = 0.05)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 预测分数等高线图\n",
    "https://scikit-learn.org/dev/auto_examples/mixture/plot_gmm_pdf.html#sphx-glr-auto-examples-mixture-plot-gmm-pdf-py\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def plot_zzz(gmm):\n",
    "    # display predicted scores by the model as a contour plot\n",
    "    x1 = np.linspace(0, 10)\n",
    "    x2 = np.linspace(0, 10)\n",
    "    x3 = np.linspace(0, 10)\n",
    "    x4 = np.linspace(0, 10)\n",
    "    x5 = np.linspace(0, 10)\n",
    "    x6 = np.linspace(0, 10)\n",
    "    x7 = np.linspace(0, 10)\n",
    "    x8 = np.linspace(0, 10)\n",
    "    X1, X2, X3, X4, X5, X6, X7, X8 = np.meshgrid(\n",
    "        x1, x2, x3, x4, x5, x6, x7, x8)\n",
    "    XX = np.array([X1.ravel(), X2.ravel(), X3.ravel(), X4.ravel(),\n",
    "                   X5.ravel(), X6.ravel(), X7.ravel(), X8.ravel()]).T\n",
    "    Z = -gmm.score_samples(XX)\n",
    "    Z = Z.reshape(X.shape)\n",
    "\n",
    "    CS = plt.contour(X, Y, Z, norm=LogNorm(vmin=1.0, vmax=1000.0),\n",
    "                     levels=np.logspace(0, 3, 10))\n",
    "    CB = plt.colorbar(CS, shrink=0.8, extend='both')\n",
    "    plt.scatter(X_train[:, 0], X_train[:, 1], .8)\n",
    "\n",
    "    plt.title('Negative log-likelihood predicted by a GMM')\n",
    "    plt.axis('tight')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "# plot_zzz(gmm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# covariances \n",
    "https://scikit-learn.org/dev/auto_examples/mixture/plot_gmm_covariances.html#sphx-glr-auto-examples-mixture-plot-gmm-covariances-py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"\\nn_components = 5\\n# Try GMMs using different types of covariances.\\nestimators = {cov_type: GaussianMixture(n_components=n_components,\\n              covariance_type=cov_type, max_iter=20, random_state=0)\\n              for cov_type in ['spherical', 'diag', 'tied', 'full']}\\n\\nfor index, (name, estimator) in enumerate(estimators.items()):\\n    # Since we have class labels for the training data, we can\\n    # initialize the GMM parameters in a supervised manner.\\n    estimator.means_init = np.array([X_train[y_train == i].mean(axis=0)\\n                                    for i in range(n_classes)])\\n\\n    # Train the other parameters using the EM algorithm.\\n    estimator.fit(X_train)\\n\\n    h = plt.subplot(2, n_estimators // 2, index + 1)\\n    make_ellipses(estimator, h)\\n\""
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'''\n",
    "n_components = 5\n",
    "# Try GMMs using different types of covariances.\n",
    "estimators = {cov_type: GaussianMixture(n_components=n_components,\n",
    "              covariance_type=cov_type, max_iter=20, random_state=0)\n",
    "              for cov_type in ['spherical', 'diag', 'tied', 'full']}\n",
    "\n",
    "for index, (name, estimator) in enumerate(estimators.items()):\n",
    "    # Since we have class labels for the training data, we can\n",
    "    # initialize the GMM parameters in a supervised manner.\n",
    "    estimator.means_init = np.array([X_train[y_train == i].mean(axis=0)\n",
    "                                    for i in range(n_classes)])\n",
    "\n",
    "    # Train the other parameters using the EM algorithm.\n",
    "    estimator.fit(X_train)\n",
    "\n",
    "    h = plt.subplot(2, n_estimators // 2, index + 1)\n",
    "    make_ellipses(estimator, h)\n",
    "'''    \n",
    "# plt.figure()\n",
    "# h = plt.subplot(1, 1, 1)\n",
    "# make_ellipses(estimators, h)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "266px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "position": {
    "height": "385px",
    "left": "765.8px",
    "right": "20px",
    "top": "110px",
    "width": "519.4px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
